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We present the results of direct numerical simulations of heavy particle transport in 
homogeneous, isotropic, fully developed turbulence, up to resolution 512 3 (R\ as 185). 
Following the trajectories of up to 120 million particles with Stokes numbers, St, in the 
range from 0.16 to 3.5 we are able to characterize in full detail the statistics of particle 
acceleration. We show that: (i) The root-mean-squared acceleration a lms sharply falls 
off from the fluid tracer value already at quite small Stokes numbers; (ii) At a given 
St the normalised acceleration a rms /(e 3 /v) 1 / 4 increases with R\ consistently with the 
trend observed for fluid tracers; (Hi) The tails of the probability density function of the 
normalised acceleration a/a rms decrease with St. Two concurrent mechanisms lead to 
the above results: preferential concentration of particles, very effective at small St, and 
filtering induced by the particle response time, that takes over at larger St. 



1. Introduction 

Small impurities like dust, droplets or bubbles suspended in an incompressible flow are 
finite-size particles whose density may differ from that of the underlying fluid, and cannot 
thus be modelled as point-like tracers. The description of their motion must account for 
inertia whence the name inertial particles. At long times particles concentrate on singular 
sets evolving with the fluid motion, leading to the apparition of a strong spatial inhomo- 
geneity dubbed preferential concentration. At the experimental level such inhomogeneities 
have been long known (see lEaton fc Fessler 19941 for a review) and utilised for flow vi- 
sualisation (e.g. exploiting bubble clustering inside vortex filaments). The statistical de- 
scription of particle concentration is at present a largely open question with many in- 
dustrial and environmental applications. We mention spray combustion in Diesel engines 

IjPost fc Abraham 2002j) or some rock et propellers (|yiI|ejiieu L jj^j^^ the for- 

mation of rain droplets in warm clouds ( |Pinsky <k Khain 1997| |Falkovic h, Fouxon fc Stepanov 2 002 
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Shaw 2003 ) or the coexistence of plankton species ^Rothschild fc Osborn 1988l|Lewis fe Pedley 2000| >. 
Inertial particles are also relevant to spore, pollen, dust or chemicals dispersion in the 
atmosphere where the diffusion by air turbulence may be even overcome by preferential 
clustering ( |Csanady 1980|lSeinfeldT986|l . 

From the experimental side, the study of particle motion in turbulence has recently 
undergone rapid progress thanks to the development of effective optical and acousti- 
cal tracking techniques IjLa Porta et al 20011 ILa Porta et al 20021 IMordant et al 20011 
|Warhaft, Gylfason, fc Ayyalasomayajula 20051). In para llel with experimental effort, the- 
oretical analysis "^aT^oro^^^^I^oyic^ fc Fouxon 200TllFalkovic¥lF Pumir 20041|Bec , Gawedzki fc fT 
Zafohflt 1 jlimojii^^ (|Boivin, Simo7rmT^SqTirres^998| 

Reade fc Collins 2000l|Zhou, Wexler fc Wang 2001|IChun et al. 2005|) are paving the way 
to a thorough understanding of inertial particle dynamics in turbulent flows. Recently, the 
presence of strong inhomogeneities characterised by fractal and multifractal properties 
have been predicted, and found in theoretical and numerical studies of stochastic laminar 
flows HBalkovsky, Falkovich fc Fouxon 200 j] |Bec, Gawedzki fc Horvai 2004| IBec 2005j) . 
in two dimensional turbulent flows jBolfeTtaTl^e'TnTo'Tr'Gamba 2004| | and in three di- 
mensional turbulent flows at moderate Reynolds numbers in the limit of vanishing inertia 
UFalkovich fc Pumir 2004) . 

Here we present a Direct Numerical Simulations (DNS) study of particles much heavier 
than the carrier fluid in high-resolution turbulent flows. In particular, we shall focus on 
the behaviour of particle acceleration at varying both Stokes and Reynolds numbers. 
For fluid tracers, it is known that trapping into vortex filaments (|La Porta et al. 20011 
Bife rale et al. 2005)) is the main source of strong acceleration events. On the other hand, 
little is known about the acceleration statistics of heavy particles in turbulent flows, 
where preferential concentration may play a crucial role. Moreover, since in most applied 
cases it is almost impossible to perform DNS of particle transport in realistic settings, it 
is very important to understand acceleration statistics for building stochastic models of 
particle motion with and without inertia (Sawfo rd fc Guest 19 911. 

The material is organised as follows. In § 2, we briefly recall the equations of motion 
of the inertial particles and summarise the DNS set up. In § 3, we present and discuss 
the main results concerning the acceleration statistics. In § 4 we discuss the acceleration 
statistics conditioned on the local structure of the flow. The final Section is devoted to 
conclusions and perspectives. 

2. Heavy particle dynamics and numerical simulations 

The equations of motion of a small, rigid, spherical particle immersed in an incompress- 
ible flow have been consistently derived from first principles by |Maxey fc Riley 1983| In 
the limiting case of particles much heavier than the surrounding fluid, these equations 
take the particularly simple form 



Here, X(t) denotes the particle trajectory, V(t) its velocity, u(x,t) is the fluid velocity. 
The Stokes response time is t s = 2p p a 2 /(9pfu) where a is the particle radius p p and pf 
are the particle and fluid density, respectively, and v is the fluid kinematical viscosity. 
The Stokes number is defined as St — Tg/r^ where t v = (V/e) 1 / 2 is the Kolmogorov 
timescale and e the average rate of energy injection. Eq. (|2.1() is valid for very dilute 
suspensions, where particle-particle interactions (collisions) and hydrodynamic coupling 
are not taken into account. 




dV 

dt 



V(t)-u(X(t),t) 



(2.1) 
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Table 1. Parameters of DNS. Microscale Reynolds number R\, root-mean-square velocity 
u rms , energy dissipation e, viscosity v, Kolmogorov lengthscale n = (/v 3 /e) 1//4 , integral scale 
L, large-eddy Eulerian turnover time Te = L/u Ims , Kolmogorov timescale t v , total integration 
time Ttot, duration of the transient regime T tr , grid spacing Ax, resolution TV 3 , number of tra- 
jectories of inertial particles for each Stokes Nt saved at frequency t v /10, number of particles 
N p per Stokes stored at frequency 10t^, total number of advected particles Ntot- Errors on all 
statistically fluctuating quantities are of the order of 10%. 



The fluid evolves according to the incompressible Navier-Stokes equations 

— +u-Vu = - + vAu + f, (2.2) 

dt Pf 

where p is the pressure field and / is the external energy source, (f ■ u) = e. 

The Navier-Stokes equations are solved on a cubic grid of size N 3 for N = 128, 256, 512 
with periodic boundary conditions. Energy is injected by keeping constant the spectral 
content of the two smallest wavenumber shells IjChen et al. 1993jl . The viscosity is chosen 
so to have a Kolmogorov lengthscale r\ » Ax where Ax is the grid spacing: this choice 
ensures a good resolution of the small-scale velocity dynamics. We use a fully dealiased 
pseudospectral algorithm with 2 nd order Adam-Bashforth time-stepping. The Reynolds 
numbers achieved are in the range R\ £ [65 : 185]. 

The equations of fluid motion are integrated until the system reaches a statistically 
steady state. Then, particles are seeded with homogeneously distributed initial positions 
and velocities equal to the local fluid velocity. Equations (|2.1(l and (12.21) are then advanced 
in parallel. A transient in particle dynamics follows, about 2 — 3 large scale eddy turn over 
time, before reaching a Lagrangian stationary statistics. It is only after this relaxation 
stage has completely elapsed that the real measurement starts. We followed 15 sets of 
inertial particles with Stokes numbers from 0.16 to 3.5. For each set, we saved the position 
and the velocity of Nt particles every dt — l/10r^ with a maximum number of recorded 
trajectories of Nt = 5 • 10 5 for the highest resolution. Along these trajectories we also 
stored the velocity of the carrier fluid. At a lower frequency ~ lOr^, we saved the positions 
and velocities of a larger number N p of particles (up to 7.5 • 10 6 per St at the highest 
resolution) together with the Eulerian velocity field. We have also followed fluid tracers 
(St = 0), that evolve according to the dynamics 

^ = «(*(<), t), (2.3) 

in order to systematically assess the importance of the phenomenon of preferential con- 
centration at varying both St and R\. 

A summary of the various physical parameters is given in table ^ 



3. Results and discussion 

In this paper we focus on the statistics of particle acceleration a(t) — From 
previous studies on fluid tracers we know that acceleration statistics is very intermit- 
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Figure 1. (a) The normalised acceleration variance a rms /(e 3 /^) 1//4 as a function of the Stokes 
number for R\ = 185 (□); R\ = 105 (o); R x — 65 (A). The inhomogeneous distribution 
of particle is quantified for the highest Reynolds in the inset, where we plot the correlation 
dimension, Z?2, as a function of St. The correlation dimension is defined as p(r) ~ r° 2 (for r <^J7 ) 
where p(r) is the probability to find two particles at distance smaller than r IBec et al. 200511 . 
(b) Comparison between the acceleration variance, a rms (□), as a function of Stokes, with the 
acceleration of the fluid tracer measured on the particle position, ((-^f-) 2 } 1 '' 2 (+)• The last curve 
(o), approaching the a rms for large St, is the one obtained from the filtered tracer trajectories, 
af mB . All data refer to Re\ — 185. 



tent and strong fluctuations are associated to trapping events within vortex filaments 
IjLa Porta et al. 2001llLa Porta et al. 2002llMordant e£ al 2001IIBiferale et al. 2005jl . How 
does inertia impacts acceleration statistics? A good starting point to gain insight on the 
effect of inertia is given by the formal solution of Eqs. (|2.1|) in the statistically stationary 
state, relating the instantaneous particle velocity to the previous history of fluid velocity 
along the particle trajectory. The expression is 

V(t) = - [ e-^- s y r 'u(X(s),s)ds (3.1) 

yielding for the acceleration 

1 r* 

a{t) = — e^^~ s ^^[u{X{t),t)-u{X{s),s)]ds . (3.2) 

J — oo 

It is instructive to analyse separately the two limiting cases of small and large Stokes 
numbers. 

At small St, i.e. r s <C t v , the fluid velocity along the trajectory evolves smoothly in 
time and the above expression for the acceleration reduces to a(t) ~ 4fU(X(t),t), i.e. 
to the derivative of fluid velocity along the inertial particle trajectory. At sufficiently 
small St this is indistinguishable from the fluid acceleration ^{X(t),t) evaluated at 
particle positions. The latter, in turn, is essentially dominated by the — Vp contribution. 
Therefore we are led to draw the following picture for the small St case: the heavy particle 
acceleration essentially coincides with the fluid acceleration; however, inertial particles arc 
not homogeneously distributed in the flow and concentrate preferentially inside regions 
with relatively small pressure gradient (low vorticity regions). As a result, the net effect 
of inertia is a drastic reduction of the root-mean-squared acceleration a rms = (a 2 ) 1 / 2 , 
due essentially to preferential concentration. Indeed, as shown in Fig.[IJi the acceleration 
variance drops off very fast already at quite small St values. In Fig. we give evidence 
that the value of a rms is very close for St < 0.4 to ((-^f-) 2 ) 1 / 2 when the average is not 
taken homogeneously in space but conditioned to be on the same spatial positions of the 
inertial particles. The agreement of the two curves supports the arguments above. Notice 
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St (a) 
<fi 2 > 
<2 4 > 



0.16 0.27 0.37 0.48 0.59 0.69 0.80 0.91 1.01 1.12 1.34 1.60 2.03 2.67 3.31 
3.09 2.07 1.80 1.63 1.50 1.39 1.31 1.24 1.17 1.12 1.06 0.97 0.88 0.75 0.61 0.51 
288 48.1 30.5 22.4 17.7 14.5 12.3 10.6 9.20 8.11 7.21 5.77 4.47 3.11 1.94 1.29 



St< b > 
<5 2 > 



St (c) 
(a 2 ) 
(a 4 ) 



0.16 0.27 0.38 0.49 0.60 0.71 0.82 0.93 1.04 1.15 1.37 1.64 2.08 2.74 3.40 

2.63 1.89 1.65 1.45 1.38 1.29 1.21 1.14 1.08 1.03 0.98 0.89 0.80 0.68 0.54 0.45 

133 32.9 21.6 16.3 13.1 10.9 9.29 8.03 7.01 6.18 5.48 4.37 3.36 2.23 1.39 0.90 

0.16 0.26 0.37 0.47 0.58 0.68 0.79 0.89 1.00 1.10 1.31 1.57 1.99 2.62 3.25 

2.02 1.59 1.40 1.28 1.19 1.11 1.05 0.99 0.94 0.89 0.85 0.77 0.70 0.59 0.47 0.39 

52.8 19.1 13.1 10.1 8.24 6.95 6.01 5.24 4.61 4.11 3.67 2.95 2.32 1.59 0.97 0.63 



Table 2. Normalised values of the second and fourth moments of the acceleration 
(a 2 ) = (a 2 )/[3(e 3 /^) 1/2 ] > <5 4 ) = i^) /[3(e 3 /v)] for (a) Rx = 185, m Rx = 105 and (c) Rx = 65. 
The statistical error on all entries are of the order of 5%. 



that at increasing Stokes the two curves start to deviate from each other, the tracer 
acceleration conditioned on the particle positions has a minimum for St w 0.5 close to 
the maximum of clustering (sec inset of Fig. ^i), eventually recovering the value of a rms of 
the unconditioned tracers for larger St. The latter effect is a clear indication that inertial 
particles explore the small scale structures of the flow more and more homogeneously by 
increasing St. In this limit a different mechanism is responsible for the reduction of the 

^rms- 

At large St, i.e. t s 3> t v , the inspection of Eq. (|3.2|1 shows that the main effect of 
inertia on particle acceleration is a low-pass filtering of fluid velocity differences, with a 
suppression of fast frequencies above t~ . In figure \\]p we also compare the acceleration 
variance with the one obtained by an artificial low-pass filtering based only on the fluid 
tracers trajectories. For each tracer trajectory, x(t), we define a new velocity, u F , filtered 
on a window-size of the same order of the Stokes time: 

1 /"* 

u F (t) = - e- < - t - s V T °u(x{s),s)ds (3.3) 

The filtered acceleration is thus given by a F = -fjjU F . Of course, in order to extract 
the effect due to filtering only we are compelled to employ fluid trajectories: (|3.3|l ap- 
plied along particle trajectories is nothing but Eq. I|3.1|l . so that the acceleration would 
coincide with the particle acceleration by definition. The root mean square fluctuation, 
a rms = ((^ mF ) 2 ) 1 ^ 2 j is tnus computed by averaging along the tracer trajectories without 
any condition on their spatial positions, i.e. homogeneously distributed in the whole 3d 
domain. The curves corresponding to a rms and to af ms become closer and closer as St 
grows larger, supporting the conjecture that preferential concentration for St > 1 be- 
comes less important. For intermediate St we expect a non trivial interplay between the 
two above mechanisms that makes very difficult to build up a model able to reproduce 
even the qualitative behaviour. 

Another interesting aspect shown in Fig. ^ is the residual dependence of the nor- 
malised particle acceleration on Reynolds number. For the case of fluid tracers it is known 
that intermittent corrections to the dimensional estimate a rms = ao(e 3 /^) 1 ^ 4 may explain 
the Reynolds dependence IjSawford et a,l. 20031 iHill 20021 iBiferale eV a,l. 2004j) . Data sug- 
gest that the fluid intermittency may be responsible of such deviations at St > as well. 
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Figure 2. (a) Acceleration pdf's for a subset of Stokes values 
(St = 0,0.16,0.37,0.58,1.01,2.03,3.31 from top to bottom) at R x = 185. The inset dis- 
plays the acceleration flatness, (a 4 )/(a 2 ) 2 , at increasing R\ from bottom to top. (b) The two 
more external curves correspond to the acceleration pdf for St = 0.16 (□) and the pdf of the 
fluid tracers acceleration measured on the same position of the inertial particles, -2^ (solid 
line). The two inner curves are the acceleration pdf at the highest Stokes, St = 3.31, (o) and 
the pdf of the filtered fluid acceleration (solid line). All curves are normalised to have unit 
variance. 



This view is supported by the fact that the curves for the three Reynolds numbers are 
almost parallel. 

A two-parameters formula for the variance of the acceleration as a function of Stokes 
number can be derived in the limit of vanishing Stokes numbers as: a% ms (St) = aj? ms (0) + 
Cexp [— (D/St) s ] (jFalkovich 20Q5J). This expression follows from the acceleration pdf of 
tracer particles under the assumptions that (i) the main effect of inertia is to reduce the 
particle concentration in regions where the accelerations is larger than i/ 1 / 2 /r a 3 / 2 ; (ii) 
the pdf tail is well reproduced by a stretched exponential shape with exponent (3 = 2/ 36. 
Although the formula fits well the data, the limitation of our data-set to only a few points 
with St <C 1 does not permit a significant benchmark of the model. 
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In table El we summarise the values that we have measured for (a 2 ) and (a 4 ) as a 
function of all Stokes and for all Reynolds numbers available. Besides the effect of inertia 
on typical particle accelerations it is also interesting to investigate the effects on the form 
of the probability distribution function a(t). As shown in Fig. Et-, the pdf's get less and 
less intermittent as St increases. In the inset of the same figure we show the flatness, 
(a 4 ) I (a 2 ) 2 , as a function of St. The abrupt decreasing for St > is even more evident 
here (notice that the y scale is logarithmic). 

In the limits of small and large St the qualitative trend of the pdf's can be captured by 
the same arguments used for a rms . In Fig. [2b we compare the pdf shape for the smallest 
Stokes number with the one obtained by using the tracer acceleration measured on the 
particle position, As one can see the two functions overlap perfectly, confirming that 
the only difference between fluid particles and inertia particles for small Stokes is due 
to preferential concentration. In the same figure we also compare for the highest Stokes, 
St = 3.31, the pdf of the particle acceleration with the one obtained from the filtered 
fluid trajectories. Now the agreement is less perfect but still fairly good, reassuring that 
this limit can be captured starting from a low-pass filter of fluid tracer velocities. It is 
worth mentioning that the pdf of tracer acceleration measured on the particle position, 
, approaches the unconditioned pdf as St increases (not shown) . This further confirms 
that preferential concentration has a minor role on the acceleration at these large Stokes 
values. 



4. Statistics of acceleration conditioned on the flow topology 

We now focus on particle acceleration statistics conditioned on the topological prop- 
erties of carrier flow at particle positions. In particular, we look at the sign of the dis- 
criminant (see e.g. |Chong, Perry fc Cantwell 1990| and lBec 2005|l : 




being &ij — diUj the strain matrix evaluated at the particle position X. Note that, 
in deriving (|4.1|l . we omitted the term proportional to Tr[<j] because of incompressibil- 
ity For A < the strain matrix has 3 real eigenvalues (strain dominated regions), for 
A > it has a real eigenvalue and 2 complex conjugate ones (rotational regions). For 
a similar study, using a different different characterisation of the flow structures, see 
also ( |Squires &l Eaton 1991| ). Note that in two-dimension the equivalent of A is the well 
known Okubo- Weiss parameter that discerns elliptic from hyperbolic regions of the flow. 

In Figs.|3^,b,c we show the acceleration pdf, P(a\A), conditioned on the sign of A at 
particle positions, for three different characteristic Stokes numbers St — 0.16,0.48, 1.34. 
In FigG3i we show the root mean squared acceleration, y/ (a 2 1 A)/3, as a function of St. 
A few results are worth to be commented. The fraction of particles in the two regions 
(iV(A ^ 0)) varies considerably as a function of the Stokes number (see inset of Fig.|3Ji), 
with a depletion of particles in the regions with some degree of rotation, which however 
becomes less effective at large St. This is similar to what is observed in the inset of 
Fig.Q., where the non- homogeneous particle distribution is characterised in terms of the 
correlation dimension (|Bec et al. 20051) . Further, despite the shape of the pdf for a given 
Stokes number does not change much as a function of the sign of A, a noticeable change 
in the squared acceleration is observed. As shown in Fig.|21i, the acceleration is higher in 
the strain dominated regions than in the ones with some degree of rotation. We remark 
that the effect of inertia is dramatic: for the smallest St the conditional acceleration is 
larger when A < while the opposite behaviour is observed for tracer (St = 0). This 
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Figure 3. Acceleration statistics conditioned on the sign of the discriminant A defined in 14. H . 
(a) Pdf of acceleration for St = 0.16 conditioned on strain regions (solid line, A < 0) and on 
rotating regions (points, A > 0) regions, respectively, (b) and (c) same as (a) for St — 0.48 and 
St — 1.34. (d) Normalised root mean square conditional acceleration on A < (empty boxes) 
and A > (full boxes) regions as a function of St. The inset displays the fraction of particles in 
the rotating regions N(A > 0) (iV(A < 0) = 1 - N(A > 0)) as a function of St. The conditional 
acceleration was computed on the data recorded at frequency IOt,, (see table 0. For St — the 
acceleration \J (o 2 |A) /3 is estimated by using the pressure gradient — Vp. 



may be the signature of the expulsion of particles out of intense vortex filaments (which 
is more effective for St <C 1) leading to an undersampling of the acceleration in the 
regions dominated by rotational motion. The same difference is also measured for higher 
moments of the conditioned acceleration (not shown). 

These results point out that the strong correlation between flow structure and particle 
preferential concentration is more effective at low Stokes numbers. At larger St the 
particle fraction A^(A ^ 0) approaches the tracer value (the response time is too large 
to maitain the correlation between particle trajectories and the local flow topology) and 
the depletion of acceleration should be ascribed to the effect of filtering, as discussed in 
the previous section (c/r. Fig.^3 and^Jj). 



5. Conclusions and perspectives 

A systematic study of the acceleration statistics of heavy particles in turbulent flows, at 
changing both Stokes and Reynolds numbers has been presented. The main conclusions 
are (i) preferential concentration plays an almost singular role at small Stokes. Indeed, 
even a quite small inertia may suffice to expel particles from those turbulent regions 
(vortex cores) where the most intermittent and strong acceleration fluctuations would 
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have been experienced; (ii) for small Stokes, a good quantitative agreement between the 
inertial particle acceleration and the conditioned fluid tracer acceleration is obtained; (iii) 
at large Stokes, the main effects is filtering of the velocity induced by the response Stokes 
times. For St > 1 , the statistical properties of fluid tracers averaged over a time window 
of the order of r s are in a quite good agreement with the inertial particle properties. 
Some important questions remain open. 

It is not clear how to build up a phenomenological model that is able to describe the 
inertial particles acceleration as a function of both Stokes and Reynolds numbers. For 
example, a naive generalisation of the multifractal description, successfully used for fluid 
tracers l|Biferale et al. 20 04). may be insufficient. In fact, it is not straightforward to 
include in such models the correlation between preferential concentration and the local 
topological properties of the carrier flow. Here such correlations have been studied in 
terms of the real or complex nature of the eigenvalues of the strain matrix at particle 
positions. We found that, more effectively at small St with respect to larger St values, 
particles preferentially concentrate in strain dominated regions and that this preferential 
concentration has a clear role in determining the acceleration fluctuations. However, this 
information does not directly bring to a model for the acceleration statistics. 

The strong fluctuations of both Kolmogorov time and Kolmogorov dissipative scale are 
certainly the most interesting aspects which distinguish the statistics of heavy particles 
in turbulence from the one measured in smooth flows. It would be then important to 
study also the statistical properties conditioned to the local Stokes number (defined in 
terms of a "local" energy dissipation, see e.g. ICollins fc Keswani 2004|l . 

Work in this direction will be reported elsewhere. 

We acknowledge useful discussions with G. Falkovich, E. Bodenschatz and Z. Warhaft. 
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